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Abstract 

The extension of the density functional theory (DFT) to include pairing correlations 
without formal violation of the particle-number conservation condition is described. 
This version of the theory can be considered as a foundation of the application of ex- 
isting DFT plus pairing approaches to atoms, molecules, ultracooled and magnetically 
trapped atomic Fermi gases, and atomic nuclei where the number of particles is exactly 
conserved. The connection with the Hartree-Fock-Bogoliubov theory is discussed. The 
method of the quasilocal reduction of the nonlocal theory is described. This quasilocal 
reduction allows to obtain equations of motion which are much more simple for the 
numerical solution than the equations corresponding to the nonlocal case. 
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1 INTRODUCTION 



Although both the theory of superconductivity and the density functional theory (DFT) have 
a long history, models which take into account pairing correlations within the framework of 
the DFT have appeared not so long ago. The initial version of the DFT developed in the 
early papers of Hohenberg, Kohn, and Sham |2I did not include the pairing correlations 
explicitly. The first generalization of the DFT in this direction was developed in Ref. [2] 
for superconductors. However, the essential feature of this theory is the nonconservation of 
the number of particles in the superconducting Fermi system. Consequently, the application 
of the DFT for superconductors to atoms, molecules, ultracooled and magnetically trapped 
atomic Fermi gases, and atomic nuclei, where the pairing correlations may be important 
but the number of particles is exactly conserved, requires an additional foundation. In fact, 
the same question emerges in connection with the DFT plus pairing approaches based, for 
example, on the theory of finite Fermi systems of Migdal (see |lj and references therein) 
or on the local density approximation (Ref. jS]). So the first goal of the present paper is 
a rigorous formulation of the extended version of the DFT taking into account the pairing 
correlations under the condition of particle-number conservation. 

The second goal is to extend the DFT to the case of a functional dependence on the 
total nonlocal single-particle density matrix (DM). Such an extension is especially important 
in applications to nuclei because it allows to introduce in a natural way the kinetic-energy 
and the spin densities dependence of the energy functional. That kind of dependence leads 
to the appearance of a radial-dependent effective mass and a spin-orbit potential which are 
essential components of the nuclear structure models. Let us note that the original version 
of the DFT can be classed as a local theory because in the papers 12] the energy functional 
only depends on the local particle density. One of the possible ways of a nonlocal extension 
of the DFT was considered in Ref. [0]. However, this method, which can be referred to as a 
straightforward extension, faces with serious difficulties related with the equations of motion 
and their physical interpretation. Another method has been recently developed in Ref. ^ to 
extend the DFT by considering an energy functional which depends on a Slater-determinant 
DM. This approach leads to the quasilocal density functional theory and allows to avoid the 
difficulties arising from the method described in jH] but, in general, the resulting DM is not 
the exact DM because only its diagonal part is the quantity which coincides with the exact 
local particle density of the interacting fermion system. 
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In the present paper we will show that the inclusion of the pairing within the frame- 
work of the extended DFT is enough to obtain the exact total nonlocal DM on the one 
hand, and to avoid the difficulties encountered in the equations of motion of Ref. [B] on 
the other hand. Although we develop our formalism in the particular case of atomic nu- 
cleus, the main results and conclusions can also be applied to any Fermi system with a 
fixed number of particles. The paper is organized as follows. In Sec. |21the extended density 
matrix (EDM) formalism is revised. In Sec. 01 the DFT is extended in order to include the 
pairing correlations. In Sec. IH the reduction to an extended quasilocal density functional 
theory is performed. The conclusions are given in the last section. In the Appendices some 
mathematical details and comments are presented. 

2 EXTENDED DENSITY MATRIX FORMALISM 

Despite the formalism of the density matrix extended to include pairing correlations in the 
ground state is well known (see, e. g., Refs. [El Ej), we shall draw some basic definitions 
which are necessary for the further analysis. Let \l/ be some arbitrary antisymmetrized 
many-fermion wave function. In the general case \E' is assumed to be normalized, but it is 
not supposed to be an eigenfunction of the particle-number operator. Thus the normal (p) 
and the anomalous (k) density matrices in the state are defined through the following 
expectation values: 

p{x,x') = {-^\a'^{x')a{x)\-^) , (1) 

k(x,x') = (^1 a(x') a(a;) 1^) , k*{x,x') = {"^la^x) a\x') \<iJ) , (2) 

where a^{x) and a{x) are creation and annihilation operators of particles in the coordinate 
representation of the usual single-particle space. In the case of atomic nuclei, symbol x = 
{r, a, q} includes the spatial coordinate r and the spin projection a variables as well as the 
index q = n,p indicating the nucleon type (neutrons and protons). 

It is important to note that for Fermi systems, which are considered here, the ground 
state (GS) is described by a wave function \I' = "^^^ with a fixed number of particles. 
Therefore, in the GS the anomalous DM vanishes: 

/tQg(x, x') = . (3) 

However, in many physical problems one can construct auxiliary quantities, which have the 
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sense of an anomalous DM but nevertheless take nonzero values in the GS even if the Eq. Q 
is fulfilled. 

Let us suppose that the DM p defined by Eq. is given for some wave function ^ 
(in particular, it may be the exact GS wave function ^^g)- If ^ is time-reversal invariant, 
we can introduce the canonical basis (CB) {^^(a;)}, where the single-particle multiindex A 
contains the sign of the spin projection s and the set of the remaining quantum numbers {c} 
(A = {c, s}), so that the following expansion is fulfilled 

c,s 

In the case of real Fermi systems, this DM is not idempotent (i. e. (? ^ p\ This is due to the 
fact that all or almost all the eigenvalues of p, defined by the equation J dx' p{x, x') (j)^ ^{x') = 
v1 (f)^ (hereinafter j dx means the space integral over r and the sum over a and q indices), 
lie in the interval < < 1. If one uses, as it was done in the theory developed by 
Gilbert in Ref. [Hj, an energy functional -E^^fp] depending only on such nonlocal and non- 
idempotent DM, in the resulting equations of motion all partially occupied natural spin 
orbitals (i. e. functions <Pcs{^) which < < 1) are eigenfunctions of the single-particle 
pseudo-Hamiltonian h^{x,x') = 6E^[p]/6p{x',x) with the same eigenvalue (see Ref. |Sj for 
details). This fact leads to difficulties in the physical interpretation and in the mathematical 
foundation of the theory developed by Gilbert. To avoid this problem, first of all we shall 
define an extended density matrix (EDM) TZ which has to be idempotent (7^^ = TZ) on the 
one hand and has to contain a given DM p as a block on the other hand. An EDM means a 
DM which is defined in a space being the doubled usual single-particle space. 

Let {ipx-riix] x)} be an arbitrary set of basis functions in this extended space spanned 
by the coordinates {x,x}, where x = ^"^^ V = ^-^e additional indices introduced for 
denoting the different components of the single-particle functions. The meaning of these 
indices will be specified in the following (in particular, the index x indicates the upper and 
lower components of the functions ipx-n{x; x) according to notations of Ref. 0). The usual 
conditions of orthonormality and completeness are supposed to be satisfied: 

/ ^^^A;r,(3;;X)^A';^'(2;;X) = Sr,,ri' Sx,X' , (5) 
Y'^\;vi^'^X)i^X;v(^''^X') = Sx,X'^i^^^') ^ (6) 

where S{x, x') = S{r — r') S^y 5q,q'- In addition we also assume that the functions ipxivi^' x) 
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satisfy the condition: 

'^x;r,ix;x) =i^l;^^{x;-x) . (7) 

From the conditions (0) and © it follows that the functions ipx-n{x', x) form a unitary 
matrix in the extended space defined previously. (Strictly speaking, they form a unitary 
operator. The use of the term matrix implies that the configuration space is discretized 
and restricted by a finite number of points.) If the condition (|7j) is also fulfilled, the Bloch- 
Messiah theorem (see Refs. jHlEIl) can be applied to this matrix. In order to reformulate 
this theorem in coordinate representation let us first introduce a complete set of orthonormal 
functions {(j)x{x)}. Notice that these functions form a unitary matrix D in the notation 
of Ref. jH] according to the rule: Djfc = (px^i^x.). Second, let us divide the set of the 
single-particle indices A into three subsets: two sets of conjugate indices p and p which 
represent "paired" states and the set of the indices h corresponding to "blocked" states, i. e.: 
{A} = {p} U {p} U {h}. Let, further, v\ and ux be real non-negative numbers which satisfy 
the following conditions: Uy^ = ^/l — , Vp = Vp , < Vp < 1 , v'^ = . According to the 
Bloch- Messiah theorem the functions ipx-^riix', x) can be represented in the following form 

^A;+(a:;x) = X]C'A'A^A';+(a;;x) , i)x--{^]x) =^Cl,^i)x'--{x]x) , (8) 

A' A' 

where C^^iy^ is a unitary matrix, and the functions tpx-ri{x', x) have the form: 
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(9) 



Let us now define the EDM TZ in terms of the arbitrary set of functions {ipx-^rii^'i x)} 
introduced above by the formula 

7^(x, X ; x\ x!) = Yl ^A;-(a^; x) rx;-{^': x') ■ (lo) 

A 

Using Eqs. © - O it can be easily shown that the following equalities are fulfilled: 

n^ = n, n^ = n, (ii) 

T^ix, X ; x', x') = Sx,x' ^(a;, x') - 7^(x', -x' ; X, -x) ■ (12) 
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Let us introduce notations for the blocks of the EDM taking into account the prop- 
erties (ITTl) - (IT^: 

7l{x, + ; x', +) = p{x, x') , TZ{x, + ; x' , — ) = k{x, x') , 

7?,(x, — ; x', +) = — x') , 7^(x, — ; x', — ) = 5{x, x') — p*(x, x') , 

or in the matrix representation: 

\-K' 1 - p- j 

As it will be clear in the following (see Eqs. 1)211) - (j^^ below), the blocks p and k play the 
role of the normal and anomalous density matrices in some (quasi-particle vacuum) state. 
Substituting Eqs. (jH)) and © into Eq. (fTUj). we obtain: 

p{x,x') = ^vl^^{x)4^l{x) , (15) 

A 

where the sum over A implies the sum over the three sets of indices p, p, and h. Thus, the 
functions (pxi^) form the CB in which p{x,x') is diagonal. So in the following we shall refer 
to Eqs. (jni) as the CB representation of the functions ip\-rf{x; x)- 

Comparing Eqs. (jH) and ()15|) we see that the expansion is a particular case of 
(fT3j) . Indeed, setting the number of blocked occupied states in (fT3j) to be even or to be equal 
to zero, one can choose vx = vx, 4>x{x) = 0a(^)? which leads to the coincidence of the right- 
hand sides of Eqs. (jlj) and (fTH|) . So, there exists a variety of sets of the functions ilJx;ri{x; x), 
which enter into the definition ()1U|) and differ by an arbitrary unitary transformation of the 
type Eq. (jSI), such that the following equalities are fulfilled: 

TZ{x, + ; x', +) = p(x, x') = p(x, x') . (16) 

Due to the fact that Eq. (0} is fulfilled for any physically meaningful DM p, actually we have 
proved that for an arbitrary nonlocal DM p(x, x'), corresponding to some interacting time- 
reversal invariant fermion system, we can construct an EDM TZ which satisfies the conditions 
(jllj) - p2)l and which is related with the DM p by the formula (fTHIl . 

It is important to note that, if the normal DM p is produced by some wave function 
according to Eq. (^, the anomalous DM k produced by the same wave function according 
to Eqs. ^ does not coincide, in general, with the quantity k defined as a block of the 
EDM by Eqs. ()13)1 . even if Eqs. (jlfij) are fulfilled. In particular, if \E' = \E'^g then, as it was 
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mentioned above, n = /t^g = 0, but for the interacting system k ^ 0. On the other hand, if 
\1/ is a quasiparticle- vacuum wave function (see below), then the equahty k, = k is fulfilled. 

In order to reproduce an arbitrarily given DM p{x, x') as a block of the EDM TZ we 
have started from an arbitrary complete set of basis functions {'ipX;Ti{x] x)} which satisfies 
the conditions (0) - ©• It is useful to do it in a different way, looking at this problem 
from another more traditional point of view. To this end, let us introduce the creation and 
annihilation operators of the quasiparticles a\ and through the equation: 

"a \ /"rfx f ^^;+'^^'+^ ^A;+(a;;-) \ | a^{x) \ ^^^^ 

With this definition the functions ipxi-qix] x) form the matrix of the Bogoliubov transforma- 
tion. Then the properties ® ^ © simply follow from the requirement of the unitarity of 
this transformation and from the fact that a\ and are a Hermitian conjugate pair of the 
operators. The inverse relation, which follows from (jH]) and (jZj), reads: 

A \^A;-(a;;+) ^A;+(a;;+) 

Let us define the quasiparticle-vacuum wave function \& by the ordinary condition: 

a^|§)=0, VA. (19) 

This definition is unique, and the wave function ^ defined by the condition (fTU|) is invariant 
under the C-transformation (jHj). So, one can find the explicit general form of \E' by using the 
CB representation The result is well known (see, e. g., Ref. jS]) and reads: 

I ^) = lliH + «I) lli^P + % 4 4) I ) ' (20) 

b p 

where | ) is the particle- vacuum wave function, cl\ = J c/x 0_;^(x) a^(x). If, further, we 
formally define the EDM TZ through Eqs. ()13p with 

p{x,x') = {^\a\x')a{x)\^) , (21) 






'4] 







k{x,x') = {^\a{x')a{x)\^) , k*{x,x') = {^\ a^x) {x') \^) , (22) 

then it is easy to show using Eqs. (|T8|l and ()19|1 that this EDM is expressed in terms of 
the functions ipx;'q{x; x) by Eq. (fTn|) . and consequently it satisfies Eqs. (|TT|l - (|T^ . In other 
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words, Eqs. dHl), (EH), (jH set a mapping of ^ onto the EDM 7^ satisfying Eqs. ^ - 
()12|) . In the following the many-to-one mappings of the many-fermion wave functions (e. g., 
\E') to the density matrices (e. g., p or 7^), defined by equations of the type (fT3|) . pT]) . and 
(j22I), will be denoted as \E' ^ p or ^ ^ 7^. Due to the fact that in the above derivation 
the matrix elements of the Bogoliubov transformation iljx;r]{x',x) c'-^'s constrained only by 
conditions © - (0), the existence of the pointed mapping \^ — ^ 7^ is true, in particular, if 
ipXynix; x) are chosen in such a way that Eqs. (jlfjj) are fulfilled for some arbitrarily given DM p 
(that is always possible as was proved above). Then we have also the mapping ^ TZ p. 
Consequently, it can be argued that for an arbitrary nonlocal DM p{x,x'), corresponding 
to some interacting time-reversal invariant fermion system, including the case for which the 
number of particles is exactly conserved, there exists a quasiparticle-vacuum wave function 
such that ^ ^ p. The explicit form of is given by Eq. ((201), in which v1 and (pxi^) are 
the eigenvalues (f^) and the eigenf unctions of the given DM p. The opposite is also 

true: any quasiparticle-vacuum wave function, which satisfies the condition — p, has the 
general exphcit form (jSO)) with ux = ux, v\ = v\, (pxi^) = 0a (^)- 

These statements can be considered as a generalization of the Lieb theorem jTJ proved 
for the local particle density and the Slater-determinant wave functions. It is remarkable that 
the inclusion of the pairing is enough not only to prove a more general statement but also 
to make the proof much more simple. Moreover, the proof of Lieb is based on the particular 
example of a Slater-determinant wave function which produces a given local particle density. 
However, in the Lieb theorem the general explicit form of such Slater-determinant wave 
function is not constructed. In contrast, in our case the general explicit form ()2()j) of the 
quasiparticle-vacuum wave function, which satisfies the condition \E' ^ p, is known. 

3 EXTENSION OF THE DFT 

Let H be the nonrelativistic exact many-body Hamiltonian of an interacting fermion system. 
Let us define an auxiliary functional which depends only on the normal nonlocal DM p: 

E[p] = inf , (23) 

'I'— >p 

where \E' are arbitrary normalized many-fermion wave functions, including the ones with a 
fixed number of particles. Following the method of Ref. j2j let us introduce an effective 



8 



many-body Hamiltonian H which generally does not coincide with H . Now we define 

£[p, k,k*]= _ mi {^\H\^) , (24) 

where \& are the quasiparticle-vacuum wave functions. Due to the existence of the mapping 
\I' — s> p proved in the previous section, the functional £[p^ k, k*] can be defined for an arbi- 
trary nonlocal DM p{x, x') corresponding to some interacting time-reversal invariant fermion 
system with a fixed number of particles and for those matrices k* which are produced 
by the quasiparticle-vacuum wave functions according to Eqs. fl22|) . Notice, however, that 
Eq. ()24j) implies that H is not a completely arbitrary operator because it is constrained by 
some mathematical conditions. First, the energy functional £ has to be well defined. This 
is not a trivial property (in spite of the existing mapping \& ^ p) because in the case of 
atomic nuclei the expectation value of the exact many-body Hamiltonian H obtained with 
the quasiparticle-vacuum wave function can diverge due to the short-range singularity of the 
bare nucleon-nucleon (NN) forces. Thus, it is assumed that the effective Hamiltonian H 
contains NN forces whose matrix elements are well defined. Second, H has to be chosen to 
ensure the minimal property of the total energy functional (see Eq. (|26|) below and Ref. [2j 
for more details). Using Eqs. and we can define, in analogy with Ref. [Zj, the 
residual correlation energy E^^: 

E^^ [p] = E[p]- mf £[p, k, k*] , (25) 

and the total energy functional £: 

£[p, k, k*] = £[p, k, k*] + E^^[p] . (26) 

The main property of the total functional £[p, k, k*] is: 

inf £[p, k, k*] = inf ( inf £[p, k, k*] + E [p] ] = inf E[p] = E , (27) 

p, K, K* p \k,k* / p 

where E^^ is the exact ground-state energy of the interacting system. If the infimum of 
the total functional ()2fij) is the minimum (usual assumption which has to be fulfilled by 
an appropriate choice of the effective Hamiltonian H), it is attained for the true nonlocal 
ground-state DM p as it follows from Eq. ()23p and from the results of Sec. El 

It is advisable to pass from the variables p, k, and k* in the functional £ to the 
components of the EDM TZ using Eqs. for k and k* , and the relation 

p{x, x') = ]^[ 5{x\ x) - nix', - ; x, -) + n{x, + ; x', +) ] (28) 



which follows from Eqs. (fT^ at p = p. Taking into account these relations, we introduce the 
energy functional S^^^ depending on the EDM: 

^ext[^]=^[A^'^*]- (29) 

It is obvious from Eqs. and that mifi£^^^^[JZ] = E^^. 

In order to establish the equations of motion of the theory, which will in the following 
be referred to as the Extended Density Matrix Functional Theory (EDMFT), let us define 
the functional 



rx;-] = S^JT^] + ^J2e, fdx 1^, 

- - dxdx' p{x,x')[6{x',x) + ^x'J^{x', x; X, x)], (30) 

X 

where Ex and p{x,x') = Pq6{x,x') are Lagrange multipliers introduced to ensure the nor- 
malization condition for ipx--{x', x) (see Eq. (jSJ) and the neutron and proton numbers con- 
servation: 

^ dr p{r, a, q; r, a, q) = Ng (31) 

which are introduced through Eq. p8|) . Applying the variational principle to the functional 
F, we obtain the following set of equations of motion (see Appendix A for details): 



^ j dx'n{x, x; x', x')i^x-r,ix'] x') =VEx1pX;v(^'i X) 



(32) 



where 

'5TZ{x\ x!\ X, x) 
These equations can also be written in the matrix form 



n{x, x\ x', x') = 2 ,^,^/^;i [ - X Sx,x' Kx, x') . (33) 



where h is the single-particle pseudo-Hamiltonian, A is the operator of the pairing field: 

_ 6S[p, k, k*] _ 6S[p, k, k*] 

' ^(^'^)- 2 , (35j 

and it is denoted: i})'^^^ = ilJx;r){x;±). Obviously, we can consider Ex > (if Ex < 0, 
the permutation iIjx;+{x; x) ^ ^A;-(a^;x) is made). Thus, the Lagrange multipliers pg and 
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E\ play the role of the chemical potential for the nucleons of the type q and the role of 
the absolute value of the quasiparticle energy respectively. From Eq. the sense of the 
indices r] and which appear in the functions ip\-ri{x; x) introduced in Sec. 121 can be easily 
understood. The index r] is the sign of the eigenvalue in Eq. (jH^ . and x denotes the upper 
and the lower components of the eigenfunctions. In most cases one of these components is 
small, and it completely vanishes for the blocked states in the CB representation © of the 
functions ipx-r){x',x)- However, in general, the functions of the CB representation, that is 
ipxyn{x; x), are not solutions of Eq. (jS2I), i- e. they are not eigenfunctions of the operator 7i 
(see comments in the Appendix B). 

As it can be seen from Eq. (jHU, the equations of motion of the EDMFT have the 
same form as the ones of the Hartree-Fock-Bogoliubov (HFB) theory. However, the EDMFT 
is not reduced to this theory. First of all, the total energy functional ()2fjj) of the EDMFT 
has a more general form as compared with the HFB energy functional. Actually, only the 
part S corresponds to the HFB functional, whereas the term E^^ has, in general, a more 
complicated and less obvious functional dependence on the DM p. But the basic difference 
between the HFB theory and the EDMFT is determined by the following: the EDMFT is 
(in principle) an exact theory in the sense that Eqs. ()27p are fulfilled. This is ensured by the 
fact that the functional E[p\, which enters in the term E^^, is defined by Eq. (j23|) through 
a set of wave functions \1/ which contains the exact wave function \l/„„ with a Gxed number 
of particles. In this sense one can say that in the EDMFT the number of particles is exactly 
conserved, despite the auxiliary quantity R, which has the sense of an anomalous DM (but 
which does not coincide with Kq^), takes nonzero values. In this context, the HFB theory can 
be considered as the phenomenological realization of the EDMFT. The relationship between 
EDMFT and HFB approaches is analogous to the relationship between the theory developed 
in Ref. [3| and the density-dependent Hartree-Fock theory, as it was discussed in more detail 
in Ref. 0. 

4 REDUCTION TO THE EXTENDED QUASILOCAL THEORY 

The theory developed in the previous sections is essentially a nonlocal one as can be seen from 
Eqs. fl34j) and (jSSI)- The exact solution of that equations is a rather complicated problem. 
However, one can noticeably simplify the task of solving these equations by reducing the 
total energy functional S[p, k, k*] to a quasilocal form following the method described in 
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Ref. [71. 

In the particular case of atomic nucleus let us introduce a set of local quantities 
consisting of the local particle n^, the kinetic-energy r,, and the spin Jq densities for neutrons 
and protons which are obtained from the nonlocal DM p as: 

n,M = ^ .'), (36) 

'^qi'^) = / dx5{x^x){Vr^r')p{x.,x), (37) 

Jqir) = dx'6{r-r')6q,q,[{cT)^,,^xVr]pix,x'). (38) 

cr 

Note that in contrast to the analogous definitions of Ref. [7 , in this case the exact DM 
p enters in Eqs. (jS^ - (jHEI) (see remark after Eq. (I2Z|)), so the quantities and Jq are 
the exact (correlated) neutron and proton kinetic-energy and spin densities respectively. 
Explicit expressions for these local quantities, which follow from Eqs. and from the 
representation of the EDM in the form are: 



riq 



^) = EEl^A;-(r",a,g;+)r, (39) 

a A 

= EEl^^A;-(^,^,9;+)r , (40) 

a A 

Jqi'^) = ^$^$^^I;-(^,^'>?;+)[('^)a',a X V]^;,._(r,a,g;+). (41) 

(t,(t' a 

We can also define the quantities which have the meaning of local anomalous densities 
for each kind of nucleons as: 

^9(^) = dx'6{r ~ r')6q^q'{ay)„f^„k{x,x') (42) 

cr 

= i^^'4^l;-{r,a',q]-){ay)„>^„ijy^_{r,a,q]+) . (43) 

<t,(t' a 

Let us now introduce three kinds of the short notations: p^^ = Up, Tn, Tp, J„, Jp}, 
X = {x„, Xp}, and x* = {x*, x*}. Using these notations, we define the following quasilocal 
energy functionals: 

s:, = inf S[p,kr^*], (44) 
S^^^Ip^^, >i, X*] = inf mf ini ^S[p, k, k*]. (45) 
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Note that the mappings p — ^ Pql, k ^ and n* x* are estabhshed according to Eqs. 
(jSni) - dSHl) and (jlD). From the definitions (glj) - (gS)), it immediately follows that the 
property (|27jl is also true for both energy functionals £^^^ and E'^^'^. Namely, we have 

p H-j'^'^'Ki^^ = p ^"i..^''''[^QL' ^' = ^GS • (46) 

In the following we shall refer to the theories based on the functionals and ()45|) as the 
extended quasilocal density functional theories EQLTl and EQLT2. 

The equations of motion for these quasilocal theories have the same matrix form (jSlJ, 
as for the nonlocal theory, but with different definitions of the operators h and A. Making 
use of Eqs. (jH^j) and (jHUj) - we obtain 



h{x, x') 



V ( [-VrT^-^Vr + f/,(r)l 5.,., - tW,{r) ■ [Vr X (T],A5{r - r') , (47) 
y\_ 2m* [r) J J 



where for the quasilocal theory EQLTl: 

^2 5^QLi <5£QLi <5£QLi 

and analogous expressions with the replacement of S'^^^ by S^^^ in the case of the quasilocal 
theory EQLT2. The difference between these two versions of the quasilocal approach consists 
in the definition of the operator A. Within the EQLTl we have 

A{x,x') = -2 ^ (49) 
dK*{x', X) 

as in the nonlocal theory, whereas in the case of the EQLT2 the operator A is purely local 
and it is obtained using Eq. 



§gQL2 

A(x, x') = -2t6{r- r') 5,,,, (a,),,, -—- . (50) 

Let us stress that the solution of the equations of motion associated with both EQLTl 
and EQLT2 enables us (at least in principle) to calculate the exact values of all the local 
densities entering in the set p^^ . The theory of Ref . allows to find the exact values of the 
local particle densities nq{r) only. This difference comes from different character of the many- 
fermion wave functions used for the building of the energy functional: a Slater-determinant 
wave function in j?! for which only the Lieb theorem was proved, and a quasiparticle- vacuum 
wave function in the present paper which enables one, as it was shown, to reproduce an 
arbitrary nonlocal DM p. 

13 



Comparing our approach with theory developed in Ref. f^, it should be noted that 
we do not introduce any external anomalous pair potential as it was done in [Sj. In our 
method the pairing emerges only as a consequence of the interaction between the fermions. 
In absence of the interaction the pairing field A vanishes, which is in agreement with the 
particle-number conservation condition. The same difference exists between the equations 
of motion of the EQLT2 case and the Bogoliubov-de Gennes equations ^2] in which the pair 
potential enters, in fact, as an external field (though created initially by the interaction). 

It should also be pointed out that the EDMFT in its nonlocal and quasilocal versions 
can be reduced to the more simple DFT plus BCS description of the pairing correlations if 
the operator 7i in Eq. ()32|) is (or is assumed to be) diagonal in the CB representation (jH)). 
However, in general, this description is not equivalent to the one obtained by solving exactly 
the EDMFT equations of motion (see comments in the Appendix B). This reduction of the 
EDMFT is analogous to the replacement of the HFB equations by coupled Hartree-Fock 
plus BCS equations how it is discussed in Refs. |Hlini- 

5 SUMMARY AND CONCLUSIONS 

We have extended the recently developed quasilocal density functional theory to include 
pairing correlations. This new approach named Extended Density Matrix Functional Theory 
is based on an extended density matrix (EDM) formalism. The EDM TZ contains as a block 
the normal density matrix p and as another block an auxiliary quantity k which has the 
sense of an anomalous density matrix, but which does not coincide with the exact one k, in 
general. The matrix k is chosen in such a way that the equality 7^2 = 7^ is fulfilled. The EDM 
which possesses this property can be easily constructed for a given density matrix p using the 
canonical basis in which p is diagonal. It has been shown that for an arbitrary density matrix 
p, corresponding to some interacting time-reversal invariant fermion system, there exist a 
qusiparticle-vacuum wave function and an EDM TZ such that the following many-to-one 
mappings take place: ^ TZ —>■ p. This statement can be considered as a generalization 
of the Lieb theorem. Using the connection between TZ and p, we have defined the total 
energy functional as an extended functional of TZ. It has been proved that its minimum 
value is equal to the exact ground state energy of the considered system. This extended 
energy functional is reduced further to a quasilocal form. Thus, in the corresponding final 
equations of motion the single-particle pseudo-Hamiltonian h and the pairing field A are 
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both (quasi) local. Although the equations of motion in our theory have the same form as 
the ones of the Hartree-Fock-Bogoliubov theory, the main difference is that the Extended 
Density Matrix Functional Theory is, in principle, an exact theory in the sense that the true 
ground stated energy can be reached for the true normal density matrix p. Finally, notice 
that the general formalism developed in this work, which introduces the pairing correlations 
in the DFT, has been discussed in the particular case of the atomic nucleus, but it can be 
easily applied to other Fermi systems with a fixed number of particles. 
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APPENDIX A 

In this Appendix the equations of motion of the EDMFT are derived. Varying the functional 
F defined by Eq. ()3()j) and taking into account Eq. (fTIH) . we obtain: 

5F = I dx dx' I ^ ^l/^'^^l^ ^) - ^ ^x.x' [XK^. x) - ^(x, x') E, ] | 

This leads to the following equations of motion: 

j dx"H{x, x; x', x') i^x;-ix'; x') = -Ex^Px;^ix; x) , (A2) 

x' 

where the operator Ti is defined by Eq. fj33p . The Eq. (A2) formally defines only one half of 
the complete set of the eigenfunctions of H. In order to define second half, let us note first 
that if the following equalities are fulfilled 

T-(-{x, x; x\ x') = -n{x', -x \ x, -x) = n*{x', x'; x, x) , (A3) 



15 



then the functions ip^.j^lx; x), defined through the functions x) using the condition 

((Tj), are also the eigenfunctions of H. In this case the complete set of the eigenfunctions of 
Ti. is divided into two equal parts with eigenvalues +Ex and —Ex- On the other hand, if the 
condition (jZj) is fulfilled, and if the set of the functions ^px■-i^'■: x) is taken as a half of the 
complete set {ipx-±{^i x)}? then Eq. (fT^ is true, and it follows from Eq. that Eqs. (A3) 
are also fulfilled. Consequently, properties (|7j) and (A3) follow from each other, and there 
exists a solution of the equations of motion which possesses both properties. So, setting 
Eq. (0) to be satisfied, we actually choose a solution which has the symmetry defined by this 
equation without imposing additional constraints in the variational procedure. In that case, 
using Eqs. ((Tj) and (A3), we come from (A2) to Eq. (j221) which represents complete set of 
the equations of motion. 

APPENDIX B 

In this Appendix the relationship between the eigenfunctions of the operators TZ and 7i is 
analyzed. First of all, notice that from Eqs. (fTUI) and (0) it follows that 

j dx'n{x, X; X', X)'^x-rj{x'; X) = Sr,,-'lpX;rji^', X) ■ (Bl) 

x' 

Thus, the functions ipx;r}{x; x) are eigenfunctions of the operator 71. Obviously, the set of the 
eigenfunctions of TZ is determined by a given TZ up to an arbitrary unitary transformation 
of the type Eq. (jH)). Consequently, the functions ipx-r]{x;x) defined by Eqs. (jHI) are also 
eigenfunctions of the operator TZ. On the other hand, from Eqs. (Bl) and ()32j) it follows 
that there exists at least one set of eigenfunctions (namely the set {ipx;-q{x',x)}) which is 
common for both operators TZ and Ti. However, from these equations it does not follow that 
any eigenfunction of TZ, and in particular ipx;-q{x] x), will be also eigenfunction of 7i. Indeed, 
while a set of eigenfunctions of TZ is determined up to an arbitrary unitary transformation 
(jSJ, this is not true for the eigenfunctions of 7i, since Eqs. and (jH^ are not covariant 
under this transformation. Consequently, if we use an arbitrarily given representation of 
the eigenfunctions of the operator TZ, in particular the CB representation Q, an additional 
nontrivial C-transformation, which does not change TZ, is generally required to diagonalize 
Ti (see Refs. |H1 E] for details). 
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